% fig2a.m ——— 复刻论文 Fig.2(a)：x 峰值分岔图 ———
clear; clc; close all;

%% —— 参数设置 —— 
beta   = 1;      % 论文中 β
xi     = 0.8;    % 论文中 ξ
mu     = 1;      % 论文中 μ
A      = 3.4;    % 外激励幅值 A
omega  = 1;      % 外激励角频率 ω

dt     = 0.01;   % 时间步长
T0     = 200;   % 丢弃前 T0 时间
Tsim   = 800;   % 采样 Tsim 时间
N0     = floor(T0/dt);
Nsim   = floor(Tsim/dt);

S0 = [0.2; 0.1; 0.2];   % 初始状态

alpha_min = 1.0;
alpha_max = 4.0;
Nalpha    = 600;
alpha_list = linspace(alpha_min, alpha_max, Nalpha);

peak_x = cell(1, Nalpha);

%% —— 并行循环提取峰值 —— 
parfor ia = 1:Nalpha
    alpha = alpha_list(ia);
    S = S0;
    x_ser = zeros(N0+Nsim,1);
    for k = 1:(N0+Nsim)
        t = (k-1)*dt;
        mu_s = A * cos(omega*t);
        % RK4 迭代一步
        k1 = f(S, alpha, beta, xi, mu, mu_s);
        k2 = f(S+0.5*dt*k1, alpha, beta, xi, mu, mu_s);
        k3 = f(S+0.5*dt*k2, alpha, beta, xi, mu, mu_s);
        k4 = f(S+    dt*k3, alpha, beta, xi, mu, mu_s);
        S = S + dt*(k1+2*k2+2*k3+k4)/6;
        x_ser(k) = S(1);
    end
    % 丢弃瞬态，仅对后续数据提峰
    x_cut = x_ser(N0+1:end);
    pks = findpeaks(x_cut);
    peak_x{ia} = pks;
end

%% —— 绘制分岔图 —— 
figure;
hold on;
for ia = 1:Nalpha
    a = alpha_list(ia);
    scatter( repmat(a, size(peak_x{ia})), peak_x{ia}, 5, 'b.' );
end
xlabel('\alpha','FontSize',12);
ylabel('x 峰值','FontSize',12);
title('Fig.2(a) 分岔图','FontSize',14);
xlim([alpha_min alpha_max]);
grid on;





